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Abstract. In this paper we present detailed models of the vertical structure (temperature and density) of passive 
irradiated circumstellar disks around T Tauri and Herbig Ae/Be stars. In contrast to earlier work, we use full 
frequency- and angle-dependent radiative transfer instead of the usual moment equations. We find that this 
improvement of the radiative transfer has strong influence on the resulting vertical structure of the disk, with 
differences in temperature as large as 70%. However, the spectral energy distribution (SED) is only mildly affected 
by this change. In fact, the SED compares reasonably well with that of improved versions of the Chiang & Goldreich 
(CG) model. This shows that the latter is a reasonable model for the SED, in spite of its simplicity. It also shows 
that from the SED alone, little can be learned about the vertical structure of a passive circumstellar disk. The 
molecular line emission from these disks is more sensitive to the vertical temperature and density structure, and 
we show as an example how the intensity and profiles of various CO lines depend on the adopted disk model. 
The models presented in this paper can also serve as the basis of theoretical studies of e.g. dust coagulation and 
settling in disks. 
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1. Introduction 

The dust continuum emission often observed from T Tauri 
stars and Herbig Ae/Be stars is widely believed to orig- 
inate from a dust/gas disk surrounding these stars (see 
e.g. Beckwith & Sargent jl996 ). These disks presumably 
have mass between 1O~^M0 and fewxlO~^MQ, and are 
thought to be the disks of dust and gas from which the 
stars were formed. At early stages of the disk's evolution, 
the disks may still accrete, and lose mass to the central 
star (Calvet, Hartmann, & Strom Calvet et al. (2000| )). 
The energy released during this process is then the main 
source of power for the dust continuum emission from the 
disk. At later stages of the disk's evolution the mass accre- 
tion rate drops, and the outer regions of the disk start to 
become dominated by stellar irradiation instead of viscous 
energy release. As the accretion rate drops even further, 
more and more of the disk becomes "passive" , and irradi- 
ation eventually will be the only source of heating for th e 
disk. As was discussed by Kenyon & Hartmann ( 1987 ), 
such a passive disk can adopt a flaring shape, i.e. a shape 
in which the ratio of the vertical surface height Hg (R) over 
R increases with R, the distance from the central star. In 
this way the disk's surface always "sees" the stellar sur- 



face, and therefore intercepts some fraction of the stellar 
radiation at all radii. 

The detailed vertical structure of such a disk can be 
computed by solving the equations of vertical pressure 
balance coupled to the equations of radiative transfer. 
Such computations have been done by e.g. D'Alessio et 



al. ( |I998| , |I999D , and Bell et al. ( |1997|J1999| ). Their mod- 
els include detailed physics, including viscous dissipation, 
heating by cosmic rays, and a study of the effects of 
self-gravity. A basic weakness of these models is their 
simplified treatment of radiative transfer, which uses the 
frequency-integrated moment equations in the Eddington 
approximation with Planck- and Rosseland-mean opac- 
ities. The reason for this simplification is that solving 
the full angle- and frequency-dependent radiative trans- 
fer equations is a challenging technical problem. It is 
known that straightforward methods for radiative transfer 
(e.g. Lambda Iteration and Monte Carlo methods) con- 
verge extremely slowly at high optical depth, and one 
is never completely sure that a true solution has been 
reached. 

In the field of stellar atmospheres these problems are 
well known, and many different sophisticated radiative 
transfer algorithms have been developed over the years. 
Such algorithms include Accelerated Lambda Iteration, 
Complete Linearization methods, and Variable Eddington 
Factor methods (for a review of radiative transfer tech- 
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niques, see e.g. Kudritzki & Hummer 199C; Hubeny 1992; 
Mihalas |1978D . 

In particular the latter method has already been 
used successfully in application to accretion disk models. 
Hubeny ( 1990 ) showed that when the Eddington factors 
and mean opacities are known, the disk temperature at 
every optical depth can be expressed by a simple ana- 
lytic formula. He presented an analytical solution for the 
case in which the Eddington approximation is adopted 
and the Planck- and Rosseland mean opacities are used. 



Malbe t & Bertout ( |1991| ) and Malbet, Lachaume & Monin 
( |200l| ) also presented disk equations based on the variable 
Eddington factor method, and present solutions to the full 
angle-dependent transfer problem. Their solutions, how- 
ever, are based on the assumption of a grey opacity, and 
therefore fall short of our aims. A full frequency-angle de- 
pendent vertical structure model has been presented by 



Sincell & Krolik (1997) for X-ray irradiated accretion disks 



around active galactic nuclei. Technically, that work comes 
close to the kind of calculations we wish to make in this 
paper, but in addition to the different physics involved, 
their calculations only solve the upper irradiated skin of 
the disk. 

It is the goal of this paper to present complete vertical 
structure models for irradiated passive flaring disks based 
on full angle-frequency-dependent radiative transfer and 
vertical hydrostatic equilibrium. We use the method of 
Variable Eddington Factors as our main radiative trans- 
fer algorithm, and we present a variant of this algorithm 
that works fast and is stable under all circumstances. An 
Accelerated Lambda Iteration (ALI) algorithm is used to 
check the results and make sure that a true solution to 
the transfer equation is obtained. 



2. The disk model 

The structure equations for a passive irradiated disk can 
be grouped in three sets. First we have the equations de- 
scribing the transfer of primary (stellar) photons as they 
move from the star radially outwards, and eventually get 
absorbed by the dust grains in the upper layers of the 
disk. The energy absorbed in these surface layers will be 
re-emitted at infrared wavelengths. Half of this radiation 
will travel upwards and escape from the disk. The other 
half will move downwards into the disk, where it will once 
again be absorbed and re-emitted. Since the disk's optical 
depth can be rather high, in particular at short wave- 
lengths, this process can repeat itself many times. In this 
way the energy diffuses all the way down to the equato- 
rial plane of the disk. This process of radiative diffusion 
can be well described by the equations for plane-parallel 
1-D vertical radiative transfer. It is here that most treat- 
ments in the literature use the moment equations with the 
Eddington approximation and mean opacities (henceforth 
MEMO method). In this paper we treat this problem in 
a fully angle-frequency dependent way. Once this prob- 
lem has been solved and the temperature stratification 



has been found, we can proceed to solve the third and 
final equation: the equation of vertical pressure balance. 

These three coupled sets of equations are solved iter- 
atively: we move from stage 1 (primary stellar radiation), 
to stage 2 (diffuse radiation field) to stage 3 (vertical pres- 
sure balance), and back to stage 1. This is repeated until 
the relative difference in the density between successive it- 
erations drops below the convergence criterion, generally 
taken to be 10"^. 

2.1. Stage 1: Extinction of primary photons 

The impinging of primary (stellar) radiation onto the sur- 
face of the disk is modeled using a "grazing angle" recipe, 
similar to that used by D'Alessio et al. (1998) and Chiang 
& Goldreich ( |1997| , henceforth CG97). Dehberately we do 
not use full 2-D/3-D ray-tracing, since this procedure can 
lead to numerical instabilities (DuUemond in prep.). The 
"grazing angle" recipe models the irradiation using ver- 
tical plane-parallel radiative transfer, with the radiation 
entering the disk under an angle P{R) with respect to the 
disk's surface. At each radius R and each vertical height 
z we evaluate the local primary stellar flux: 



Fu{R.z) 



exp(-r,(i?,z)//3(i?)) 



(1) 



where is the stellar luminosity and the vertical optical 
depth Tiy(i?, z) is defined as 



,{R,z) 



p{R, z)K^dz , 



(2) 



with p{R, z) the dust density (g cm~'^) and Hi, the dust 
absorption opacity (cm^ g~^)- We ignore scattering, and 
we ignore the higher order geometrical effects which start 
to play a role when z/R > 1. The grazing angle P{R) is 
defined as /3(i?) = O.AR^/R + Rd[HjR)/dR, where Hs 
is the surface height of the disk and i?, the stellar ra- 
dius (see CG97). We define Hg as the height above the 
midplane where 63% (i.e., 1 — exp(— 1)) of the integrated 
stellar radiation has been absorbed. Although this defini- 
tion is somewhat arbitrary, it is not critical for the results 
presented in this paper. 

It is convenient to express I3{R) as 



/3(i?) = 0.4^+e(i?)§ 

where ^ is the flaring index defined as: 
d\g{H,/R) 



MR) 



(3) 



(4) 



This flaring index is a number of order ^ ~ 2/7 (CG97). 
Its value is computed self-consistently during the iteration 
procedure. Usually the iteration is started with a guess 
for ^ (for instance ^ = 2/7), and updated after each few 
iteration steps. In order to avoid numerical instabilities, ^ 
is always evaluated two radial gridpoints away from the 
point where it is used (see appendix of Chiang et al. 2001). 
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It is important to note that a proper self-consistent 
computation of the flaring index ^ is crucial if one wants 
to achieve energy conservation. Assuming it to be some 
fixed value is guaranteed to result in disks that emit more 
(or less) radiation than they receive. 

2.2. Stage 2: Vertical radiative transfer 

Once the function F^{R, z) is known, one can compute 
the amount of absorbed primary radiation per volume el- 
ement: 



qiR,z) 



p{R, z)k^F^{R, z)dv 



(5) 



This energy will then be re-emitted as infrared radiation, 
half of which will diffuse towards lower z into the disk inte- 
rior. The transfer of this re-emitted IR radiation through 
the disk can be approximated as a 1-D slab geometry 
transfer problem of high optical depth. 

The intensity of this diffuse infrared radiation field is 
denoted by I^ ^, and obeys the following radiative transfer 
equation: 



dh 



dz 



(6) 



where Ki, is the opacity per unit of matter, p is the material 
density and Bi,{T) is the Planck function. This differential 
equation has to be integrated along z for all v and /i, in 
order to obtain the intensity function I^^^[z). The dust 
temperature T is determined by assuming thermodynamic 
equilibrium with the radiation field: 

/ pK^B^{T)dv ^ I pK^J^dv + — , (7) 
Jo Jo 47r 



where the mean intensity is defined as 



(8) 



Equations form a set of coupled integro- 

differential equations in the /i, z space. A straightfor- 
ward way to solve them would be to use the method of 
Lambda Iteration (LI): first evaluate Ia^u{z) (Eq. ||), then 
Ji,{z) (Eq. |) and finally T{z) (Eq. gj, and iterate this 
procedure until convergence is reached. However, it is well 
known that this method converges very slowly at high op- 
tical depths (see e.g. Rybicki & Hummer 1991). 

Instead, we use the method of variable Eddington 
factors (VEF, see appendix ^ for a description of the 
method). This method estimates the shape of the radi- 
ation field by using a moment equation, but eventually 
reaches a solution that solves the full frequency- and angle- 
dependent transfer problem. It converges generally within 
~10 iterations, and is therefore much faster and more re- 
liable than the ALI method. 



2.3. Stage 3: vertical hydrostatic equilibrium 

Once the temperature at all locations in the disk is known, 
one can integrate the pressure balance equation to find the 
density structure. The vertical pressure balance equation 



dP{R, z) 
dz 



= -p{R,z) 



GM 



(9) 



Since P — kpT/prriu (with p — 2.3 for a H2, He mixture), 
and the derivatives of T are known, this equation can be 
cast into the form of an integral in p. This integration 
is started at z = 0, taking p{R, 0) at first as an arbitrary 
value. One can then compute the resulting vertical surface 
density: 



S = 2 



p{R, z)dz , 



(10) 



and then renormalize the density profile to the actual sur- 
face density. This new density structure is then compared 
to the density structure of the previous iteration. If the 
convergence criterion is not met, the iteration proceeds by 
going back to Stage 1. This iterative procedure has proven 
to be stable and to converge quickly, typically after about 
five to eight iterations. 

3. Resulting vertical structure 

As our example case we take a star with Teff — 3000K, 
i?* = 2.0 i?o and — 0.5 M©. Our disk has an inner 
radius of i?in = 3-R*, an outer radius of i?out ~ 300 AU 
and a gas -f dust surface density S = Sq (i?/AU)^^ with Sq 
the surface density at 1 AU. We assume that the gas and 
the dust are always mixed in the mass ratio 100 : 1, and 
we take for the dust opacity the opacity of astronomical 
silicate of Draine & Lee ( 1984 ) for grains of 0.1 pm size. 

In Fig. |l| the vertical structure of the disk at a dis- 
tance of lAU from the star is shown. For this calculation 
we took So = lO'^g/cm^, which corresponds to a visual 
optical depth through the disk ry = 2.3 x 10"*. To better 
illustrate the results, we compare then to those obtained 
for the same parameters using the more approximated 
MEMO method of solution. The left panel shows the ver- 
tical temperature profile. One sees that the full transfer 
model has a much smoother structure than the MEMO 
solution, and that the midplane temperature is signifi- 
cantly lower, at least for the high optical depth case shown 
here. The main effect of the lower midplane temperature 
is a slightly smaller pressure scale height at the equator 
{Hp/R = 0.022 for the VEF models versus Hp/R = 0.028 
for the MEMO one, as expected from the factor 1.57 differ- 
ence in the midplane temperature between the two mod- 
els) , so that the disk is a bit more compressed toward the 
equatorial plane. 

The right panel of Fig. shows the run of the density 
as function of z/R. One can see how the two solutions 
differ only by a small factor near the midplane, since the 
pressure scale heights are not very different. At higher z/R 
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Fig. 1. The vertical structure of a passive irradiated flaring disk surrounding a star with Tcs = 3000K, i?* = 2Rq 
and A/» = 0.5Mq. The distance from the star at which the structure is shown is 1 AU. The gas+dust surface density 
at this radius is E(IAU) = lO'^g/cm^. The dashed line shows the structure as computed using the moment method 
with mean opacities (the MEMO equations). The solid hue shows the structure computed using full angle- frequency 
dependent radiative transfer (using the VEF method). Left panel shows the temperature, while the right panel shows 
the density. The tickmarks on the curves indicate the location of the defined disk surface, i.e. the z/R above which 
63% of the direct stellar light has been absorbed. 




Fig. 2. The radiation field in the disk, as computed using the moment equations (MEMO method: dashed line) and 
the full radiative transfer (VEF method: solid line). The model and the location of the slice (1 AU) are the same as in 
Fig.^. Left panel the frequency-integrated mean intensity J and Eddington fiux H. Right panel the frequency-averaged 
Eddington factor /. In the MEMO method this factor is fixed by definition to 1/3, while in the VEF method it is 
computed sclf-consistcntly. Note that the slight deviation of / from 1/3 close to the equatorial plane is due to the 
discretization of the radiation field in /x, so that the integration of /i^ over /i is not exactly 1/3. 



the differences become much larger (an order of magnitude 
at z/R ~ 0.075 — 0.1). This can be understood as a slight 
vertical shift in z of the steep density profile. 

Fig. ^ illustrates some technical aspects of the ra- 
diation transfer solutions, which clarify the results just 
shown. Firstly, one can see (right panel) how the 
Eddington factor f{z) drops below 1/3 in the upper layers 
of the disk. This is because in these optically thin layers, 
most of the intensity is more or less parallel to the disk 
instead of pointing out of the disk. This is an effect that 



is opposite to what happens for instance in stellar winds, 
where most of the radiation eventually gets beamed out- 
wards, and the value of / becomes close to 1. 

It is also clear from Fig. ^ (left panel) that as a result 
of the use of the Rosseland mean opacity in the MEMO 
equations, the mean intensity in the MEMO method re- 
mains virtually constant, whereas in reality the mean in- 
tensity drops strongly towards the equatorial plane. This 
results in a non-isothermal disk structure, as opposed to 
what would be expected on the basis of the usual diffusion 
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approximation considerations. According to these consid- 
erations one would expect to have a perfectly isothermal 
disk interior since there are no sources of heat in the disk. 
Indeed, if a grey opacity is used in the full transfer model' 
as a test, the disk interior becomes isothermal as expected. 
But as soon as non-grey dust opacities are included, the 
isothermality is broken. The reason is that some radia- 
tion is able to leak out of the disk at long wavelengths, 
where the opacity is much smaller than at the wavelength 
of the bulk of the radiation. The slight positive tempera- 
ture gradient causes a downward diffusive flux at shorter 
wavelengths that exactly cancels this energy leakage, so 
that the total flux remains zero deep within the disk. 

In Fig.^ the vertical temperature structure is shown 
for different values of the radius. Near the star, where 
the disk optical depth is very large, the equatorial plane 
temperature is usually below that of the MEMO model, 
whereas at large radii, where the disk becomes optically 
thin to its own diffuse radiation, the MEMO models pre- 
dict lower temperature values than what we find. The lat- 
ter effect is related to the fact that the MEMO method 
treats the diffuse radiation field in a frequency-integrated 
way, and has therefore no information on how the energy 
is distributed over frequency. In other words: the MEMO 
method uses mean opacities based on the local dust tem- 
perature, whereas the local radiation temperature may be 
different in the optically thin case. 

In Fig.^ it can also be seen that at small radii the 
model-predicted temperature for large z / R is exactly con- 
stant. In reality we expect a small decrease of T for 
increasing z/R since the distance to the star increases 
as V-R^ -\- z"^ . But since higher order geometrical effects 
are ignored in our model this is not taken into account. 
Fortunately the gas-|-dust density at large z/R is very 
small, and this effect has no practical influence on the 
resulting spectra. 

Finally, we show in Fig.^ the radial dependence of the 
equatorial plane temperature (left panel) and of the pres- 
sure scale and surface height of the disk for the VEF and 
MEMO solutions (right panel). As was already known 
from FigJ^, the midplane temperature is lower than in 
the MEMO approximation at most radial intervals, but 
becomes higher at the outer edge. However, at the very 
inner edge it is again the same in both approaches. This 
is because at these radii the vertical optical depth of the 
disk becomes so large that even at millimeter wavelengths 
there is no possibility for the gas to cool, and the disk 
becomes isothermal. Fig.^, right panel, shows the radial 
profile of the pressure scale height, which follows approx- 
imately that of the midplane temperature, as discussed, 
and of the surface height. It is interesting to note how the 
MEMO approximation and the full transfer model give al- 
most equal values of Hs- This is a result of the higher tem- 
peratures at higher elevations above the midplane. This 
can be seen from the fact that the upper part of the tem- 
perature rise, which is where the direct stellar radiation is 
absorbed, is located at about the same z in both models. 



150 
^ 100 
I- 50 

88 
^ 60 
^ 40 
^ 20 

58 
40 

^ 20 
10 

38 
^ 20 

I- 10 






R = 3.5 AU 



R = 17.2 AU 



R = 84.2 AU 



R = 300.0 AU 



0.2 



0.4 

z/R 



0.6 



0.8 



Fig. 3. The vertical structure at different radii, as com- 
puted using the moment equations (dashed line) and using 
full angle-frequency dependent radiative transfer using the 
VEF method (solid line). On the horizontal axis the di- 
mensionless vertical height. The model parameters are the 
same as in Fig.^ 



With these differences between the simple (MEMO) 
and the complex (VEF) method, one may wonder whether 
a compromise between the two methods may give re- 
sults that are close enough to the real one. Even with 
the efficiency of the VEF method, solving the full angle- 
dependent radiative transfer problem is a time-consuming 
process. It may therefore not be very practical for e.g. the 
computation of the evolution and/or the dynamics of the 
disk, where the radiative transfer problem needs to be 
solved thousands of times. In Fig.^ we show the same 
slice as in Fig.^. Here we added the computed temperature 
profiles for two alternative methods: one in which the full 
//-dependence is solved, but with Planck- and Rosseland 
mean opacity replacing the kj and Kjy, and one in which 
the Eddington approximation is used, but with the full 
frequency-dependent opacities. 

One sees that the former does not improve much on the 
MEMO method. But the latter in fact gets very close to 
the 'right' answer. This shows that most of the problems 
with the MEMO method originate from the wrong treat- 
ment of mean opacities, while the use of the Eddington 
approximation iself {fy = 1/3 and H{z = Zynax) = 
J{z = 2;max)/V3) docs not cause major problems. We con- 
clude that an approximate method, in which frequency- 
dependent radiative transfer is done in the Eddington ap- 
proximation, is in fact a relatively accurate method. This 
method avoids the CPU-time consuming full frequency- 
angle dependent radiative transfer, and is therefore faster 
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Fig. 4. The temperature at the equatorial plane and the surface height of the disk, for the same model parameters 
as in Fig.|l|. 
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Fig. 5. Same figure as Figj^, but now compared to two ad- 
ditional methods. The dotted line is for the method using 
full /i-dependence, but the usual mean opacities (Planck 
and Rosseland) . The dot-dashed line is for the method us- 
ing the full frequency-dependent opacities (and therefore 
frequency-dependent transfer), but using the Eddington 
approximation for the angular dependence. 



than the full VEF method. Nevertheless, in this paper 
we adopt the complete VEF method, since it is still fast 
enough for our purpose, and it is exact, at least in the 1-D 
approximation made here. 

3.1. Effect on the spectral energy distribution 

The new vertical structure models have a spectral energy 
distribution (SED) that is not very different from the ones 
predicted by the vertical structures from the MEMO ap- 
proximation. In fact, they do not even differ much from 
the even simpler semi-analytic treatment of the Chiang & 
Goldreich model, if certain modifications are made to the 



latter (see Chiang et al. |200l| and DuUemond, Dominik 



Natta |2001| , henceforth DDNOl). This may be rather sur- 
prising, since the disk structure is so drastically different, 
with equatorial temperatures differing up to 70%. The ex- 
planation of this lies in the fact that the total emitted 
flux of the disk is prescribed by energy balance. At each 
radius it is expected that the disk emits both an optically 
thin component and an optically thick one. Each compo- 
nent carries half of the emitted flux, which equals the total 
absorbed flux. This means for instance that the optically 
thick component (at far-IR wavelengths) must have an ef- 
fective temperature that is the same in all models. In the 
full transfer models the effective temperature is therefore 
the same as the effective temperature of the Chiang & 
Goldreich model, even though the spectrum may not be a 
perfect blackbody. 

There are however some differences that are worth not- 
ing. The full transfer models predict higher fluxes in the 
mid and far-infrared than MEMO models. The difference 
is of the order of 30% in the wavelength range between 50 
and 500 /im for a T Tauri star, and between 50 and 100 
/im for a hotter Herbig Ae star. However, the shape of the 
SED is not very different, even if the full transfer mod- 
els fall slightly more rapidly towards the far infrared. The 
SEDs of the CG97/DDN01 models are somewhat flatter 
in the same range of wavelengths, since they are com- 
posed of two equally strong components (see the case of 
the HAe star, where the two "bumps" are clearly visible). 
This is a direct consequence of the discrete two-layered 
structure of the CG97/DDN01 model. In the fuU transfer 
model the surface layer has a continuous temperature gra- 
dient, which smoothly matches to the photosphere. This 
produces the smoother SEDs at far-IR wavelengths. The 
bumps in the SEDs of the CG97/DDN01 model are there- 
fore an artifact of the adopted two-layers simplification. 

The region aroung the 10 micron feature are virtually 
not affected by the full treatment of radiative transfer. 
The emission in this region is dominated in all models 
by emission from the optically thin surface layers, whose 
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Fig. 6. The spectral energy distribution computed 
using the MEMO method (dashed hne), using the 
CG97/DDN01 model (dotted line) and using full 
angle-frequency dependent radiative transfer with the 
VEFmethod (sohd line). For all models we use a face-on 
viewing angle {i = 0). Upper panel: SED for a T Tauri 
star with T^g = 3000K, i?, = 2Rq and AU = 0.5Mq. 
Lower panel: SED for a Herbig Ae star with Tcff — 9520K, 



R^, 



2.5i?f7, and = 2.4Mf^. For both cases the disk 



has S = 103(i?/AU)-^ Note that in the latter, we did not 
make a special treatment for the inner edge of the disk, 
as is required (DDNOl), hence the absence of the near-IR 
bump. 



properties are quite independent of what happens deep 
below in the disk, and is therefore much less affected by 
the different methods of solving the radiative transfer. 

3.2. Effects on molecular line intensities 

Knowledge of the correct vertical structure of disks (tem- 
perature and density) is of great importance for predicting 
the abundances of different molecular species and the in- 
tensity of their lines. In a separate paper we discuss the 
differences between various models in detail, in the context 
of molecular line observations from T Tauri and Herbig Ae 
stars, for which the SED is sufficiently well known (Van 
Zadelhoff & Dullemond in prep). Here we confine ourselves 
to a demonstration of the effects on the intensity of the 
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Fig. 7. The temperature of the disk at a radius of 220 AU, 
as a function of vertical Ay into the disk. Here Ay — 
corresponds to z = cxd. The tickmarks show te t — 2/3 
locations for the 3-2 line of CO, ^^CO and C^^O respec- 
tively. 



lines of CO and its two main isotopomers ^'^CO and C^^O, 
assuming that gas and dust temperatures are the same. 
The CO molecule can be used as a temperature tracer, 
and is therefore suited to probe the differences between 
the two models. 

A first insight of how line intensities depend on the disk 
model can be obtained from Fig. |^, which shows how the 
temperature varies as a function of the gas column density 
(for our template T Tauri star model and i?=220 AU). If 
we assume for simplicity that the line optical depth is pro- 
portional to the gas column density, and that the bright- 
ness temperature in the line is roughly equal to the gas 
temperature where the optical depth in the line is t ^ 2/3, 
one can immediately see that different models will predict 
the same brightness temperature for lines that are very op- 
tically thick, but that MEMO models will underestimate 
significantly the brightness temperature of more optically 
thin lines, which form more deeply in the disk where the 
predicted temperatures differ. 

For a true comparison of models to the observations a 
number of effects have to be taken into account, includ- 
ing NLTE effects, size of the telescope beam, and chemical 
abundances of the molecular species. The NLTE level pop- 
ulations have been calculate using a 2D Monte Carlo code 
(Hogerheijde & van der Tak 2000). The reason for calcu- 
lating the full NLTE level populations instead of LTE is 
due to the low densities in the upper layer of the disks 
(see van Zadelhoff et al. 2001). From these populations, 
the emitted line profiles were constructed, convolving the 
disk with a beam of 4.0", the apparent size of the source 
on the sky at a distance of 150 pc. In all cases the disk is 
assumed to be seen under an inclination of 60°. The abun- 
dance ratios with respect to H2 are assumed to be con- 
stant ([CO]/[H2]=8.xlO-^) in the entire disk. The abun- 
dances of the two isotopomers follow the isotopic ratios of 
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Table 1. The integrated line emission for the CO, ^^CO 
and C^^O transitions. All values have been derived after 
a convolution with a beam of 4.00 ". 
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Fig. 8. The molecular line emission for the 3-2 (left) and 
6-5 (right) transitions for the molecules CO, ^^CO and 
C^^O emitted by a disk at 150 pc. The disk is assumed to 
be seen under an inclination of 60°, with a beam of 4.0". 
The solid line represents the emission using the full con- 
tinuum transfer, the dashed line shows the results using 
the moments equations. 



[12C]/[13C]=60 and [i2C]/[i«C]=500, similar to the solar 
neighborhood. 

The molecular line emission results are shown in Fig. 
^ and summarized in Table |l|, in which the line emission 
integrated over velocity is given. The three isotoponiers, 
due to their differences in abundance become optically 
thick at different heights in the disk. CO reaches the t=1 
plane close to the disk surface, where temperature and 
density are very similar in the two solutions. The line pro- 
files shows only small changes between models with differ- 
ences up to 5% in the integrated line intensities. The ^'^CO 
becomes optically thick in the intermediate layer showing 
slightly larger differences between the two models (up to 
15%). C-'^^O has a optical depth of a few, and its lines form 
in the region of the disk where the VEF and MEMO solu- 
tions differ more. The differences between the two models 
are ^ 30%. The CO isotopomers can be used to measure 
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the vertical temperature profile in disks (e.g. Guilloteau & 
Dutrey 1994), and in this context the differences between 
different radiation transfer methods become significant. 

For each isotope, the difference between the more op- 
tically thick 3-2 and the less optically thick 6-5 follows the 
same pattern, as seen in Fig. However, the difference in 
the ratio of integrated intensities between models is small 
when compared with the typical observational errors (10- 
20 %). 

These calculations were performed assuming constant 
abundances throughout the disk. In reality this will not 
be the case due to dissociation of molecules in the upper 
layers by stellar and interstellar radiation and freeze-out 
of molecules on the grain surface at T < 20K (Aikawa et 
al. 2002). This latter process would enhance the difference 
between the two models as more freeze-out is expected 
from the colder MEMO model. 



4. Discussion 

The model described in this paper is relatively elemen- 
tary in the sense that all kinds of detailed microphysics 
are ignored. But within the limited physics with which 
it is defined, and the geometrical simplifications made, it 
is a reasonably accurate solution. The 1-D vertical radia- 
tive transfer and the vertical pressure balance are solved 
without any approximations other than the discretization 
itself. In this respect it is a step forwards compared to 
the models of passive disks published in the literature so 
far. The models worked out in this paper can be down- 
loaded from a website[], which also features a couple of 
one-annulus test setups for the radiative transfer problem 
in these disks. 

The advantage of ignoring all the complex and de- 
tailed microphysics is that we have a well-defined disk 
model with very few parameters. And all the features 
of this model can be explained in terms of two pieces of 
physics only: radiative transfer and hydrostatics. This sim- 
ple model can then serve as a basic model to which more 
physics can be added later, such as the inclusion of dust 
settling and coagulation, the inclusion of dust scattering, 
a treatment of internal heat processes such as viscous dis- 
sipation, etc. A study of the effect of dust scattering on 



|http : //www.mpa-garching .mpg . de/PUBLICATIOMS/DATA/ 
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the disk's structure and the outcoming spectrum is under 
way (DuUemond & Natta in prep.). 

There are however a few uncertainties concerning the 
global structure of the disk, even within the simple prob- 
lem definition used in this paper. These have to do with 
the possibility of self-shadowing. The present calculations 
start from the Ansatz that the disk has a flaring shape, so 
that the star's radiation can illuminate the disk at every 
radius. It then finds a solution that is consistent with this 
Ansatz. But in reality perhaps part of the disk might be 
non-flaring and reside in the shadow of the inner regions 
of the disk. This could be the result of the history of for- 
mation of the disk, or it may have to do with the outer 
parts of the disk becoming too optically thin to sustain 
a positive derivative of Hg/R. But a self-shadowed region 
in the disk could also be the end-product of an instability 
(Dullemond |2000| ; Chiang |2000| ). In these cases the Ansatz 
of positive flaring index is not confirmed, and the solutions 
might be different from the smooth flaring disk solutions 
presented here. In order to investigate the possibility of 
such alternative disk solutions, one must include full 2D 
radiative transfer and full 2D hydrostatics into the model 
calculations. Due to the enormous complexity of this prob- 
lem, in particular at high optical depths, we do not venture 
into this direction at present. 

The models presented in this paper describe the mid- 
dle and outer parts of a passive flaring circumstellar disk. 
The very inner parts have a different structure than the 
one described here. Very close to the star the dust has 
evaporated and one is left with a pure gas disk. If the gas 
is optically thin, the inner rim of the dusty part of the 
disk is directly exposed to the radiation of the star. It will 
therefore be much hotter than the rest of the disk, that is 
only irradiated on the surface. This inner rim will therefore 
puff up and consistute a new component in the spectral 
energy distribution of the disk (Natta et al. 2001 ). A phys- 
ical description of this inner rim, and the effects it has on 
the complete disk structure is described by Dullemond, 
Dominik & Natta ( ^OOll) . But a full 2D model for this in- 
ner rim again requires 2D radiative transfer, and is beyond 
the scope of this paper. 

5. Conclusion 

We have presented in this paper the results of calculations 
of the vertical structure of irradiated circumstellar disks 
in hydrostatic equilibrium. The main difference with pre- 
vious models (D'Alessio et al., Malbet et al. etc.) is the 
full angle- and frequency-dependent treatment of radiative 
transfer. It turns out that the frequency-dependent treat- 
ment of the problem is of crucial importance for obtaining 
the correct vertical structure (temperature and density) of 
the disk. Over most of the disk the midplane temperature 
is significantly lower than the value predicted by a simpler 
radiative transfer method that uses Planck and Rosseland 
mean opacities. The correct treatment of the angular de- 
pendence is of lesser importance. In applications where 
computation time is critical, the use of the Eddington ap- 



proximation is reasonable, as long as the full frequency 
dependence of the radiation field is taken into account. 

If the disk model is used as input for further calcula- 
tions, the need to implement a reliable radiation transfer 
method is obvious. This is, for example, the case when 
computing the expected intensity and profile of molecular 
lines, as shown in §4 for the lines of the CO isotopomers, 
which form at different depth in the disk. As more real- 
istic chemistry and line transfer models are used, the dif- 
ferences between different disk models is likely to become 
even more relevant (see van Zadelhoff and Dullemond, in 
prep.), and the VEF models should be adopted to con- 
strain physical and chemical parameters in circumstellar 
disk environments. 

The SED is also affected by the change in vertical 
structure, but the differences are not very large. In general 
the SED compares reasonably well with the results from 
both the MEMO approach and the CG97/DDN01 model. 
This means that the model of Chiang & Goldre ich (1997 ), 
with improvements described by Chiang et al. (2001) and 
Dullemond, Dominik & Natta ( 200l| ), is a fairly robust 
model to compute the SED of T Tauri stars and Herbig 
Ae/Be stars. 
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Appendix A: Radiative transfer with Variable 
Eddington Factors 

The idea at the basis of the variable Eddington factors 
method ca n be t raced back to the sixties (see e.g. Mihalas 
& Mihalas 1984 and references therein). For accretion disk 
theory the equations for this method were presented by 
Hubeny ( |l990[ ) and Malbet & Bertout (|l99l| ). Here we 
present a slight variation of this approach. The advantage 
of our method is that it can be used up to any optical 
depth even in cases where the opacity is non-grey and the 
internal energy dissipation (and therefore the net flux) 
within the disk is zero. 

Define the mean intensity J^, the Eddington flux H^, 
and the second moment of radiation K^, as follows: 



J. = - 



I^,.„^id^JL 



-1 
1 



(A.l) 



(A.2) 



(A.3) 



By multiplying the transfer equation (|6|) by resp. , and 
integrating over d^, one obtains the fc-th moment equation 
of radiative transfer. The first two moment equations are: 

^ = pn^BAT) - J,) (A.4) 
az 



dz 



(A.5) 



One can write the second moment K,j as a dimensionless 
factor times the mean intensity J^: 



where / is the Eddington factor. For isotropic radiation 
the Eddington factor is equal to / = 1/3. For purely 
beamed radiation {I^^u = for p. ^ 1) this factor is 
/ = 1. If one makes an assumption for this value, then 



Eqs.(A.4, A.5) form a closed set of equations, which can 
be solved subject to boundary conditions. Often the as- 
sumption is made that the radiation field will be more 
or less isotropic, and therefore that f^, = 1/3. This is the 
Eddington approximation, and stands at the basis of many 
approximate transfer codes. 

But if one had some way of computing (rather than 
guessing) the Eddington factor, then the moment equa- 
tions are not an approximation anymore, and in fact will 
be exactly equivalent to the full transfer equation. One 
way of doing so is to iteratively switch between the mo- 
ment equations and the real transfer equation. Start with 
a guess for fi,{z) (tak e it for instance 1/3), and solve 
the moment equations (A.4, A.5) to find the temperature. 
Then integrate the formal transfer equation (|^) using the 
current temperature. Using the resulting radiation field 
I^_u{z) one can then compute the Eddington factor fu{z). 
Once this is done, one starts all over again, using the newly 
computed Eddington factor. This process is repeated until 
a converged solution for T(z) is reached. Since in this pro- 
cedure the fv{z) is computed on the fly, one can be sure 
that this solution is a real solution, and not an approxi- 
mate one. It is important to note that if this solution is 
inserted into the full set of transfer equations (^,0,||) , then 
one will see that these are solved as well. In other words: 
the moment equations were used to find a solution to the 
real transfer equations. 



The advantage of the moment equations (A.4, A.5) is 
that they can be solved directly, using a two-point bound- 
ary value approach. In fact, the frequency- integrated mo- 
ment equations are even simpler to so lve directly. The 
frequency-integrated version of Eq.( A.4) is: 



^-^=p(n,iTf-T^~n,j), 
where Kp{T) is the Planck mean opacity and 



H 



H,,dv . 



(A.7) 



(A. 



Using the concept of energy conservation, one can replace 
the right-hand-side of Eq.(A.7) with the source term: 



dH_ 

dz 



_q_ 

in 



(A.9) 



This equation can be trivially integrated from the equator 
(starting with 77 = 0) up to z = Zmax- At that point we 
can compute the value of the frequency integrated mean 
intensity J: 



where 



(A.6) 



J(z ^ Zmax) = H{z = Zmax)/^', 



J — Jydv , 
'o 



(A.IO) 



(A.ll) 
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and ijj is the ratio of H/J as computed from the fuU ra- 
diative transfer. Using the frequency-integrated version of 
Eq.(O), 



djfJ) 
dz 



-P 



(A.12) 



(where / is the J- mean of /^), and starting from the value 
of J{z = Zniax) computed above, one can integrate back 
towards the equator to find J{z). The temperature then 
follows from 



n 



KjJ 



_<1_ 

An 



(A.13) 



where kj is the J^-mean of the opacity computed from 
the full transfer. 

By iterating on the above procedure, one can quickly 
find solutions to the transfer equation. An even quicker 
convergence can be reached by applying a linear conver- 
gence amplifier like Ng's algorithm (Ng 1974). 

Typically one needs about 70 gridpoints in z, 40 points 
in fj, and, dependent on the kind of opacity table and the 
wavelength range, about 60 points in v. It is important to 
make sure that there are enough /i gridpoints near fj. — 0, 
because in plan-parallel geometry most of the radiation is 
at small values of ii. We use a logaritmic grid in /i spanned 
between /zq — 0.01 and 1. 



It is important to note why, in Eq.(A.12), we did not 
use a flux-weighted mean opacity kh (or Rosseland opac- 
ity) multiplied by the frequency-integrated Eddington flux 
H, as was suggested for instance by Malbet & Bertout 
( |l99l| ). The reason is that in a disk without internal 
heat production, the total flux H is zero. The frequency- 
dependent flux H^, on the other hand, is non-zero, consist- 
ing of downward diffusive flux at short wavelengths and an 
equal amount of upward flux at long wavelengths. For non- 
grey opacities one therefore has a virtually zero H close to 
the equator, but a very non-zero K^H^dh' at the same 
location. An evaluation of kh = /g°° K^Hydv/ H^dv 
would therefore yield a diverging number, and the algo- 
rithm becomes unstable, whereas the equation in the form 
of Eq.( A.12 ) yields a stable algorithm. 

The method described in this section has proven to 
work well, and reaches a solution to an accuracy of 10~^ 
in temperature within 11 iterations and to 10~® within 22 
iterations, independent of the optical depth. On a Pentium 
III with 850 MHz clock frequency a typical radiative trans- 
fer problem is solved in about 5 seconds. 



ALI code. For those cases in which the ALI in fact con- 
verges to a satisfiable degree, we find that the solution are 
indeed very close to each other. 

To further convince ourselves of the validity of the 
solutions provided by the VEF method, we applied the 
VEF method to a problem with grey opacities. For such 
a problem it is known that the MEMO method should 
work reasonably well. And as expected, the VEF method 
is consistent with the MEMO results. 

A.2. Convergence 

In order to demonstrate the advantage of the VEF method 
over methods like LI and ALI, we show here the conver- 
gence history for a test problem. Consider a circumstellar 
disk which has S — lO'^g/cm^ in dust mass at 1 AU from 
a central star of Toff = 3000K and i?* = 2.0 i?0. The ver- 
tical pressure scale height of the disk is assumed to be 
Hp = 0.028 AU, and the density profile is assumed to be 
Gaussian. The flaring index equals 0.2, and the opacity ta- 
ble used is for astronomical silicate (Draine & Lee 1984) 
In Fig 



A.l 



the convergence history for the temperature 
is plotted using three different methods: LI with Ng, ALI 
with Ng and the VEF method with Ng. The convergence 
history is only shown for the first 100 iterations, within 
which only the VEF method converges. After about 400 
iterations, the ALI4-Ng method in fact also converges to 
within 1% of the solution, but the LI-|-Ng method does 
not converge even after 1000 iterations. 

It should be noted that the method of Accelerated 
Lambda Iteration has not been widely used for dust con- 
tinuum transfer. It has been developed mainly for line 
transfer, for which it has proven to be reasonably effec- 
tive. We adapted the algorithm for use with dust contin- 
uum transfer, and chose as our approximate operator the 
diagonal of the full lambda operator. 



A.l. Testing 

In order to verify the reliability of the code, a few tests are 
performed. First of all we do internal consistency checks. 
Since the VEF method computes the J and H from both 
the full transfer as well as the moment equations, one can 
a-posteriori check whether the two are the same. For all 
our models this turns out to be indeed the case. Then we 
compare the solution to the results from an independent 
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Fig. A.l. Comparison between different angle- and frequency-dependent radiative transfer methods for the stage 2 
transfer problem. Plotted is the temperature at each iteration step. Left plot is for Lambda Iteration (LI), middle 
plot is for Accelerated Lambda Iteration (ALI), right plot is for Variable Eddington Factor (VEF) method. All three 
methods use Ng acceleration, which is responsible for the 'jumps' in the convergence. 



